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Harmonic inversion techniques have been shown to be a 
powerful tool for the semiclassical quantization and analysis 
of quantum spectra of both classically integrable and chaotic 
dynamical systems. Various computational procedures have 
been proposed for this purpose. Our aim is to find out which 
method is numerically most efficient. To this end, we sum- 
marize and discuss the different techniques and compare their 
accuracies by way of two example systems. 

PACS numbers: 05.45. -a, 03.65.Sq, 02.70.-c 



I. INTRODUCTION 

Semiclassical trace formulae relate the spectra of quan- 
tum systems to the periodic orbits of the pertinent clas- 
sical systems |ij . They yield expansions of the quantum 
response function of the form 



9(E) = E 



dk 



E-E k 



-^po^ p 



(i) 



Here, Ek are the energy eigenvalues of the quantum sys- 
tem, dk their multiplicities, 5 po is the action of a clas- 
sical periodic orbit, A po an amplitude that can be cal- 
culated from classical mechanics (including phase infor- 
mation given by the Maslov index), and the sum on the 
right-hand side extends over all periodic orbits and usu- 
ally diverges for real energies E. Thus, the quantal infor- 
mation cannot be extracted directly from the semiclassi- 
cal expansion. 

One particular and widely applicable method to over- 
come the convergence problems of the periodic orbit sum 
is semiclassical quantization by harmonic inversion [^|,^) . 
By a Fourier transform of Eq. (|I|) the problem of semiclas- 
sical quantization can be recast as a harmonic inversion 
problem, viz. the extraction of the frequencies Uk = Ek/h 
and amplitudes dk from a given time signal 



C(t) = E dfce " 



(2) 



Especially intriguing, and important, are systems pos- 
sessing a classical scaling property, i.e., the classical dy- 
namics does not depend on an external scaling parame- 
ter w and the action S'po = ws po of periodic orbits varies 
linearly with w, with s po the scaled action. This is not 
a severe restriction since it covers a variety of systems, 
such as systems with homogeneous potentials, billiard 
systems, or the hydrogen atom in static external fields. 



For scaling systems the semiclassical signal which has to 
be harmonically inverted has the special form of a sum 
of S functions with peaks at positions given by the scaled 
actions s po of the periodic orbits, 



C(s) 



■^po^ ^po) 



(3) 



The frequencies of the signal (j|) are the semiclassical ap- 
proximations to the quantum eigenvalues u>k of the scal- 
ing parameter. By the same token, harmonic inversion 
of signals with the functional form (||) also plays a key 
role in the high resolution analysis of the density of states 
6{E) = J2 n — E n ) of quantum spectra, in an effort 
to extract information about the underlying classical dy- 
namics |3|-||]- 

In practical applications, the signal @ is always known 
in a finite range < s < 5 max only. The signal length 
•Smax is often fixed or at least hard to increase, e.g., for 
periodic orbit quantization of classically chaotic systems 
where the number of periodic orbits proliferates expo- 
nentially with the signal length. To obtain the optimum 
results from the harmonic inversion procedure, it is essen- 
tial to choose an algorithm which allows one to extract 
the most precise estimates for the spectral parameters 
{uJk, dk} from the signal of a given length S max . 

Various computational procedures have been proposed 
for the harmonic inversion of signals of the type (|3|) . How- 
ever, a systematic study of the relative merits and de- 
merits of the methods and a quantitative study of their 
efficiencies is still lacking. To remedy this situation, we 
summarize and discuss different techniques of harmonic 
inversion and compare their accuracies in the application 
to two simple albeit typical example systems for which 
exact trace formulae are known. The aim is to pin down 
the numerically most efficient method for harmonic in- 
version. 



II. HARMONIC INVERSION OF DELTA 
FUNCTION SIGNALS 

Due to the finite signal length S max , the analysis of the 
signal by conventional Fourier transform is limited by the 
"uncertainty principle" , which states that in a signal of 
finite length S max two frequencies can only be resolved 
if their difference is larger than 27r/5 max . The uncer- 
tainty principle can be circumvented by the application 
of high-resolution methods (|^] which extract a discrete 
set of frequencies and amplitudes without imposing any 
a priori restrictions on the frequencies tOk- However, an 
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uncertainty remains in the weaker form of the "informa- 
tional uncertainty principle" 0, which states that the 
signal length S max required to resolve the frequencies is 
given by 



> 



4irg(Lu) 



(4) 



in terms of the local average density of frequencies g(u>). 
Since this relation involves the average instead of the 
minimum spacing between frequencies, the signals can 
usually be significantly shorter than required by the 
Fourier transform. 



) is a nontrivial task 
y, the number of fre- 



Harmonic inversion of the signal ( 
for the following two reasons. First 
quencies contained in the signal is usually large, which 
implies that simple methods for solving the harmonic in- 
version problem may be numerically unstable. Secondly, 
the special functional form of a signal as a sum of 6 
functions does not fulfill the prerequisites of several es- 
tablished algorithms that C(s) should be known on an 
equidistant grid, s = nr, with n — 0, 1, 2, . . . We 
briefly review the four numerical methods which so far 
have been successfully applied to the harmonic inversion 
of signals C(s) given as a sum of 8 functions. 

Method 1: Discrete FD A powerful tool for the 
harmonic inversion of signals given on an equidistant grid 
is the filter-diagonalization (FD) method of Wall and 
Neuhauser || in the version of Mandelshtam and Taylor 
. The basic idea is to introduce a set of basis functions 
localized in frequency space and to recast the harmonic 
inversion problem as a generalized eigenvalue problem. 
For a suitable choice of the frequency window the subset 
of frequencies located in that window can be obtained 
by solving a generalized eigenvalue equation with small 
matrices. 

To evaluate the signal @ on an equidistant grid, the 8 
functions must be given an artificial width a which spans 
several grid points, i.e., a > r. This regularization can 
be achieved by convoluting the signal with a narrow, e.g., 
Gaussian function. In this form the FD method has been 
applied in Rcfs. 0,0] as a tool for semiclassical periodic 
orbit quantization. 

The convolution of the signal C(s) results in a damp- 
ing of the amplitudes dk — > d^ = dkexp(—w^a 2 /2). 
The width a of the Gaussian function should be chosen 
sufficiently small to avoid an overly strong damping, e.g., 
by setting a < |w m ax| _1 where ui max is the largest fre- 
quency of interest. To properly sample each Gaussian 
a dense grid with sufficiently small step size (t « er/3) 
is required. Therefore, the convoluted signal to be pro- 
cessed by FD usually consists of a large number of data 
points, in particular when high frequency regions of the 
signal are to be analyzed. The numerical treatment of 
such large data sets may suffer from rounding errors and 
loss of accuracy. 

Method 2: 8 -function FD The artificial smooth- 
ing of the signal can be circumvented when following a 
different approach suggested by Gremaud and Delande 



H . In contrast to Ref. Q where the signal is assumed to 
be known on an equidistant grid, Gremaud and Delande 
start from a continuous-time formulation of the FD algo- 
rithm close to the original ansatz in || . Integrals involv- 
ing the 8 function signal (^) can then easily be calculated 
and expressed in terms of periodic orbit sums. 

Method 3: Discrete DSD An alternative method 
to FD for high-resolution signal processing is decimated 
signal diagonalization (DSD), which was introduced by 
Belkic et al. || . DSD is a two-step process for harmonic 
inversion. In the first step a low-resolution frequency 
filter is applied by subjecting the signal to a discrete 
Fourier transform, selecting the Fourier components in 
the frequency interval of interest, and applying an inverse 
Fourier transform to them. The resulting band-limited 
signal contains only a small number of frequencies in the 
chosen interval and can therefore be further processed, in 
the second step, by conventional high-resolution methods 
such as, e.g., Linear Prediction or Pade Approximants 
P,^0[. DSD effectively uses, in this processing stage, 
the operational part of the discrete version of FD [Q, 
which constructs small matrices of a generalized eigen- 
value problem directly from digitized points of the band- 
limited decimated signal. The DSD technique is designed 
for signals given on an equidistant grid but can be ap- 
plied to the 8 function signal (||) after convolution with 
a narrow, e.g., Gaussian function in the same way as ex- 
plained above (see Method 1). 

The DSD method of Ref. |8| is easy to implement as 
it basically resorts to standard algorithms for discrete 
Fourier transform and matrix diagonalization. However, 
the application of the low-resolution Fourier filter in the 
first step of the method implicitly assumes periodicity 
of the signal (with period equal to the signal length) , in 
which case the DSD filter is exact. In general, of course, 
this condition is not met, so that only approximate filter- 
ing can be achieved. Therefore, DSD must be expected 
to be less accurate than FD (Method 1) for frequencies 
close to the borders of the window, or when very short 
frequency windows are chosen (see the discussion in Sec. 

Method 4: 8 -function DSD The DSD technique 

(Method 3) can be modified for a more direct applica- 
tion to the 8 function signal (|3|) without the necessity of 
convoluting the signal with a narrow, e.g., Gaussian func- 
tion. Because the signal C(s) is given as a periodic orbit 
sum of (5 functions the 'filtering' step can be performed 
analytically by replacing the discrete Fourier transform of 
Method 3 with the continuous form of the Fourier trans- 
form and expressing the integrals in terms of periodic 
orbit sums. This technique was proposed in Ref. | pl| . 
The application of the analytical filter for a rectangu- 
lar frequency window [wo — Aw,wo + Aw] results in a 
band-limited (bl) signal Jll| 



Cbi(s) = XMp° 



po 



sin [(s - s po )Aw] ^ 

?r(s - S po ) 



(5) 
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where the S functions are basically replaced with sine 
functions [sine x = (sin x) /x] . The band-limited signal 
(0) can be discretized with a small number of points 
and further processed, in the second step, by conven- 
tional harmonic inversion methods as described above 
(see Method 3). 

In practice, the band-limited signal can only be evalu- 
ated approximately because the signal is only known up 
to a finite length. Since the sine functions decay slowly 
at infinity, peaks well beyond the end of the known signal 
may have an influence on the band-limited signal points. 
Omitting contributions from the (unknown) peaks be- 
yond the limit of the given signal amounts to the as- 
sumption that the signal be zero outside the given range. 
Note that this filter differs from the low-resolution filter 
of Method 3 where the signal is implicitly assumed to be 
periodic. 

In summary, the four methods can be classified accord- 
ing to whether they are discrete-time algorithms (Meth- 
ods 1 and 3), which require a regularization of 6 function 
signals to be discretized, or continuous-time algorithms 
adapted to 6 function signals (Methods 2 and 4). Alter- 
natively, they can be classified into filter-diagonalization 
(FD) methods (Methods 1 and 2) and decimated sig- 
nal diagonalization (DSD) methods (Methods 3 and 4) 
where the low-resolution 'filtering' and high-resolution 
signal processing is performed in two separate steps. 



III. NUMERICAL EXAMPLES AND 
DISCUSSION 

To quantitatively assess the relative performances of 
the different algorithms, we present a comparison of the 
numerical accuracy achieved by all of these methods for 
two simple but archetypal examples, viz. firstly, the high- 
resolution analysis of the spectrum of the harmonic oscil- 
lator and, secondly, the search for the zeros of Riemann's 
zeta function as a mathematical model for periodic orbit 
quantization in chaotic systems. We choose these sys- 
tems because they possess exact trace formulae, so that 
the decomposition of the semiclassical signal in a sum of 
exponentials is known to be exact. 



A. Harmonic oscillator 

The one-dimensional harmonic oscillator (with fku = 

1) has energy eigenvalues E n = n + \, n = 0, 1, 2, Its 

density of states can be written as an exact trace formula 
HI: 



kE 



(6) 



n=0 



k— — c 



The right-hand side of Eq. (^|) can be interpreted as a pe- 
riodic orbit sum [in analogy to Eq. (Q)] where S% = 2nkE 



is the action of the (k times traversed) periodic orbit and 
dk = (— l) k is the periodic orbit amplitude. [The inter- 
pretation of the k = Thomas-Fermi term is special, see 
Ref. II for more details.] The high-resolution analysis 
of the quantum spectrum g(E) — J2^=o d(E — E n ) should 
yield equally spaced real frequencies 0Jk = 27rfc and am- 
plitudes dk = (— f) fc of equal magnitude. Thus, this sim- 
ple application of harmonic inversion to the extraction of 
classical periodic orbit parameters from a quantum spec- 
trum || 5| is ideally suited to compare the efficiencies 
of the different methods for the harmonic inversion of 6 
function signals. 

Since the signal is periodic with period AE = 1, 
an integer signal length would render the low-resolution 
Fourier filter of Method 3 exact. To avoid this atypical 
situation, we choose signal lengths as rational multiples 
of 7r. According to the informational uncertainty princi- 
ple (|4|) a signal length of £ max > 2 should suffice to re- 
solve the frequencies. Typically, Eq. (Q) slightly underes- 
timates the required signal length. We therefore present 
results calculated with a signal of length E max = 7r, which 
means that only three energy levels contribute to the sig- 
nal. To assess the accuracy of the results, we use the 
absolute values of the imaginary parts of the calculated 
frequencies and amplitudes as error indicators. If the 
analysis of the signal were exact, all imaginary parts 
should vanish. Therefore, an inspection of the sizes of 
the imaginary parts allows us to check the accuracy of 
the calculation. We note that this sort of accuracy test 
can be applied to all bound systems. If the exact frequen- 
cies are known, as is the case in our example systems, the 
real parts can also be compared. The errors of the real 
and imaginary parts typically are of the same order of 
magnitude and exhibit, at least qualitatively, the same 
behavior. 

Results for the harmonic inversion of the quantum 
spectrum g(E) obtained with the four methods intro- 
duced in Sec. || are presented in Figs. [I] and || for the 
imaginary parts of the frequencies u>k and amplitudes dk , 
respectively. For frequencies which appear to be miss- 
ing, imaginary parts of zero have been obtained by the 
pertinent method. From figure parts (a) to (f) the width 
Auj of the frequency filter is increased. For the applica- 
tion of Methods 1 and 3 the signal has been discretized 
with step width r = 0.002 after convolution of the sig- 
nal with a Gaussian function with width a = 0.006. In 
all cases it can be seen that the precision achieved de- 
creases to the boundaries of the frequency window. The 
reason is that none of the filtering methods is exact and 
can neither completely remove the influence of frequen- 
cies outside the window of interest nor exactly preserve 
the contributions of frequencies inside the window. For 
narrow windows, the FD methods I and 2 outperform the 
DSD algorithms 3 and 4, for wide windows the situation 
is reversed. The frequencies obtained by Methods I and 
2 are equally precise for small windows, whereas for wide 
windows Method 2 gains superiority and even competes 
with the DSD methods. In general, the distance from the 
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FIG. 1. Imaginary parts (absolute values) of the frequen- 
cies u>k extracted from a harmonic oscillator signal of length 
-E max = 71" • Symbols x, +, and denote to Methods 1 to 
4, respectively. Windows are [10 — Aw, 10 + Aw] with Aw = 
(a) 4.5, (b) 5.5, (c) 6.5, (d) 7.5, (e) 8.5, (f) 9.5. 

window boundaries where a method acquires its full pre- 
cision is smaller for the FD than for the DSD methods. 
Calculations were carried out with double precision. For 
the widest window shown, frequencies have practically 
been obtained to machine precision. 

For all methods, the precision of the amplitudes in Fig. 
^ is somewhat lower than that of the frequencies in Fig. |l| 
because the amplitudes are calculated from the eigenvec- 
tors of a generalized eigenvalue problem, which in general 
are less accurate than the eigenvalues. In particular, the 
difference in precision between the frequencies and am- 
plitudes is considerably larger for Method 1 than for any 
other method, so that even for small windows the ampli- 
tudes obtained by this technique are the least accurate 
(see the x symbols in Fig. |^). 



B. Zeros of Riemann's zeta function 

It is a peculiar feature of the harmonic oscillator sig- 
nal that the density of frequencies is constant, i.e., the 
precision of frequencies obtained from a signal of a given 
length is the same throughout the whole frequency do- 
main. However, in typical systems the density of states 
grows rapidly with the frequency, which means that a 
longer signal is required to extract higher frequencies. As 
an example of a system exhibiting this typical behavior, 
we discuss the Riemann zeta function which has served 
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FIG. 2. Same as Fig. [l] but for the imaginary parts of the 
amplitudes d^. 

as a mathematical model for periodic orbit quantization 
P,P^|. It is well known that, if the zeros of Q(z) on the 
critical line Re z — | are written as z = h — iw, the 
density of zeros on the critical line can be expressed as 



1 CO 

>h =—EE 

p m= 1 



hip 

i/2 



cos(wmlnp) , 



(7) 



where p runs over all prime numbers. Eq. (Q) is for- 
mally identical to a semiclassical trace formula with 
Spm = wmlnp corresponding to classical actions and 
Ap m = {\np)/p m / 2 corresponding to classical ampli- 
tudes. With this interpretation, the Riemann zeta func- 
tion can be used as a mathematical model for chaotic 
dynamical systems, and the Riemann zeros are obtained 
by harmonic inversion of the 5 function signal B 



(8) 



p m—l 



Unlike typical semiclassical trace formulae, Eq. (0) is ex- 
act. As the zeta function has only simple zeros, the am- 
plitudes dk obtained from the harmonic inversion of the 
signal (||) must all be equal to 1. 

In Fig. |^ we present our numerical results obtained by 
application of Methods 1 to 4 to extract the (numerically 
complex valued) Riemann zeros with Re Wk < 100. Ide- 
ally, all values Wk should be real. Therefore, the absolute 
values of the imaginary parts of the Wk can again serve 
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FIG. 3. Imaginary parts (absolute values) of locations Wk 
of zeros of the Riemann zeta function in the frequency window 
[10,100]. Symbols X, +, and denote to Methods 1 to 
4, respectively. Signal lengths are Smax = (a) 4.5, (b) 5.0, (c) 
5.5, (d) 6.0, (e) 6.5, (f) 7.0. 

as a measure for the accuracy of the harmonic inversion 
process. For the application of Methods 1 and 3 the 
signal has been discretized with step width r = 0.002 
after convolution of the signal with a Gaussian function 
with width a = 0.006. 

It is no problem to construct the signal (||) for the 
Riemann zeros up to quite large maximum values S ma , x 
because only prime numbers are used as input. How- 
ever, the periodic orbit quantization of physical systems 
usually requires a numerical periodic orbit search which 
becomes more and more expensive for longer orbits, es- 
pecially in chaotic systems, where the number of orbits 
proliferates exponentially with increasing signal length. 
Therefore, in practical periodic orbit quantizations the 
given signal length is often rather short. In Fig. |^ we 
present the results for the accuracy of the methods for 
harmonic inversion for various signal lengths, increasing 
from S max = 4.5 in Fig. |a to S max = 7.0 in Fig. |f. The 
frequency window w G [10, 100] is kept fixed. 

For a fixed signal length, the zeros up to a certain crit- 
ical value can be obtained to a constant precision. Above 
the critical frequency, the precision decreases rapidly due 
to the higher density of states. As was to be expected, for 
all methods the critical frequency increases with growing 
signal length, which means that frequencies in regions 
of higher spectral density can be resolved. Roughly, the 
critical frequency is determined by the informational un- 
certainty principle (^). In fact, it is slightly higher for 



FIG. 4. Same as Fig. |3| but for the imaginary parts of the 
multiplicities dk- 

the FD methods 1 and 2 than for the DSD methods 3 and 
4. As before, the maximum accuracy below the critical 
frequency is obtained by the DSD methods. However, 
above the critical frequency the precision yielded by the 
FD methods is higher. 

The lowest zero of the zeta function is located at 
w = 14.1347, not far above the lower boundary of the fre- 
quency window at w — 10. For the first zeros a decrease 
in accuracy due to the proximity of the boundary can be 
seen. Evidently, the influence of the boundary diminishes 
with increasing signal length. Again, it is considerably 
more pronounced for the DSD than for the FD methods. 
For the latter, it can only be seen in the shortest signals. 
With any of the four methods, the boundary effects on 
the lowest zeros can be removed if the lower boundary of 
the frequency window is decreased to w = 0. 

Fig. [I] presents results similar to those shown in Fig. 0, 
but for the imaginary parts of the multiplicities . The 
accuracy of the results achieved with the different meth- 
ods resemble those obtained for the imaginary parts of 
the frequencies Wk, with the exception of Method 1 (see 
the x symbols) which provides amplitudes with signifi- 
cantly lower precision. 

IV. CONCLUSION 

In this paper we have quantitatively determined the 
accuracies of four different algorithms for the high- 
resolution harmonic inversion of i5 function signals, by 
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applying all algorithms to two, physically motivated, ex- 
ample signals. For sufficiently long signals and broad 
frequency windows the four methods provide excellent 
results of very high accuracy, in the case of the exam- 
ples selected even close to machine precision. However, 
when either the width of the frequency filter or the sig- 
nal length is considerably reduced, the accuracy of the 
results obtained by the four methods can vary by several 
orders of magnitude. 

Our calculations show that no general clear-cut answer 
to the question "Which method is best in all physical sit- 
uations?" is possible. In practice, the window width can 
be regarded as a free parameter, i.e., it can usually be 
chosen sufficiently large to achieve good results before 
increasing computational effort or numerical instabilities 
become a restriction. The signal length, on the contrary, 
is often fixed or at least hard to increase. In such a 
case the choice of the algorithm for harmonic inversion 
of the signal will be essential to achieve the optimum re- 
sults. When the signal length is quite at the limit for 
convergence of the frequencies and amplitudes, the filter- 
diagonalization (FD) methods 1 and 2 provide superior 
accuracy compared to the decimated signal diagonaliza- 
tion (DSD) methods 3 and 4. For signals given as the 
sum of 5 functions Method 2 will often prove to be the 
method of choice. 

We conclude by noting that harmonic inversion tech- 
niques can be generalized so as to cope with the analysis 
also of multidimensional signals, with important applica- 
tions in other areas of physics |14| . The knowledge gained 
from the comparison of methods for one-dimensional har- 
monic inversion in this paper should also serve as a useful 
guide in future developments and applications of accu- 
rate and efficient algorithms for multidimensional high- 
resolution signal processing. 



[5] B. Gremaud and D. Delande, Phys. Rev. A 61, 32504 
(2000). 

[6] M. R. Wall and D. Neuhauser, J. Chem. Phys. 102, 8011 
(1995). 

[7] V. A. Mandelshtam and H. S. Taylor, J. Chem. Phys. 
107, 6756 (1997); 109, 4128 (1998) (erratum). 

[8] Dz. Belkic, P. A. Dando, J. Main, and H. S. Taylor, J. 
Chem. Phys. 113, 6542 (2000). 

[9] S. Marple, Jr., Digital Spectral Analysis with Applica- 
tions, Prentice-Hall, Englewood Cliffs, NJ, 1987. 
[10] W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. 
Flannery, Numerical Recipes, 2nd ed., (CUP, Cambridge, 
1992). 

[11] J. Main, P. A. Dando, Dz. Belkic, and H. S. Taylor, J. 
Phys. A 33, 1247 (2000). 

[12] M. Brack and R. K. Bhaduri, Semiclassical Physics 
(Addison-Wesley, Reading, MA, 1997), chapter 3. 

[13] M. V. Berry, Riemann's zeta function: A model for quan- 
tum chaos? in Quantum Chaos and Statistical Nuclear 
Physics, edited by T. H. Seligman and H. Nishioka (Lec- 
ture Notes in Physics 263) (Springer, Berlin, 1986), p. 1 

[14] V. A. Mandelshtam, J. Magn. Reson. 144, 343 (2000); 
Prog. Nucl. Magn. Reson. Spec. 38, 159 (2001). 



ACKNOWLEDGMENTS 



We thank P. A. Dando and H. S. Taylor for fruitful dis- 
cussions. This work was supported by the Deutsche For- 
schungsgemcinschaft and Deutscher Akademischcr Aus- 
tauschdienst. 



[1] M. C. Gutzwiller, Chaos in Classical and Quantum Me- 
chanics (Springer, New York, 1990). 

[2] J. Main, V. A. Mandelshtam, G. Wunner, and H. S. Tay- 
lor, Nonlinearity 11, 1015 (1998). 

[3] For review on the use of harmonic inversion techniques 
in semiclassical physics see J. Main, Phys. Rep. 316, 233 
(1999). 

[4] J. Main, V. A. Mandelshtam, and H. S. Taylor, Phys. 
Rev. Lett. 78, 4351 (1997). 



G 



